%[rp cp]=calcSnake(filename,getNewRoi) - Uses a snake to determine a region in an image
%This function applies a snake to an image using inital contour points to determine
%a region of interest.  The function returns the x and y image coordinates of the final 
%contour points.
%
%A description of the algorithm can be found at
% http://jasonheidemann.com/activecontours.aspx and is also described at 
% http://code.google.com/p/python-opencv/wiki/snakeAlgorithms.
%
%INPUT ARGUMENT
%filename - The filename of the image file to load.  It is also the path + filename 
%           prefix of the text files containing the inital contour points to use for the
%           snake alogorithm.
%getNewRoi - A logical value that is true if the user wishes to select the inital snake 
%           contour points. If you are running in Octave this must be false because Octave 
%           does not include roipoly
%
%OUTPUT ARGUMENT
%cp - A vector containing the x coordinates in the image of the final snakes contour points.  the first value will
%       be the same as the last value.   
%rp - A vector containing the y coordinates in the image of the final snakes contour points.  The first vlaue will
%       be the same as the last value.
%
%Author: Jason Heidemann - His paper is at (http://jasonheidemann.com/activecontours.aspx)
%Modified by: Paul Otto  - I just wrote createGaussianFilter, filterSeperable, and calcGradientMag since
%                           Jason could not find his functions:).
function [rp cp]=calcSnake(filename,getNewRoi)

if getNewRoi == true
    %This will not run under Octave
    [mask,cp,rp] = roipoly(img);
    save(strcat(filename,'_cp.txt'),'cp','-ascii','-tabs');
    save(strcat(filename,'_rp.txt'),'rp','-ascii','-tabs');
else
    cp = load(strcat(filename,'_cp.txt'));
    rp = load(strcat(filename,'_rp.txt'));
end

sigma = 4;
img = imread(filename);

if isrgb(img)
    img = rgb2gray(img);
end

img = im2double(img);

cp = round(cp);
rp = round(rp);
cp = round(cp);
rp = round(rp);

[H]=createGaussianFilter([sigma sigma],[11 11]);
smoothedImg=filterSeperable(H,img);

gradMag = calcGradientMag(smoothedImg);

N = size(cp,1);

%show the image with the contour points
imshow(img,[]);
hold on;
plot(cp,rp,'r');
plot(cp,rp,'bx');
hold off;
pause(0.25);
stats = [];

% Varying Sigma Mod
% allow a large sigma to direct our points to the highest gradient points
% then decrement sigma so the these our points become fixated on the actual
% highest point. So sigma is 5 for 1st iteration, 4 for 2nd, etc until
% sigma is 1
% This method seemed effective on the circle, however on the swan more
% improvement could be, possibly, if we made non-continuity penalty a lower value

for j = 1:20
    
    d = [rp-circshift(rp,1) cp-circshift(cp,1)];
    d = sqrt((sum(d.^2,2))); %L2 distance between points on the snake. 
    stats = [stats;  mean(d(2:N)) var(d(2:N))];
    for n = 1:N-1
        [newX,newY] = MinimizeE(rp,cp,n,gradMag,7,d);
        rp(n) = newX;
        cp(n) = newY;
        if n == 1
            rp(N) = newX;
            cp(N) = newY;
        end
    end
    
    figure(1)
    imshow(img,[],'notruesize');
    hold on;
    plot(cp,rp,'r');
    plot(cp,rp,'bx');
    hold off;
    
    figure(2)
    imshow(gradMag,[],'notruesize');
    hold on;
    plot(cp,rp,'r');
    plot(cp,rp,'bx');
    hold off;
    
    
    pause(0.05);
    % update sigma, and our gradient
    if sigma > 1
        sigma = sigma - 1;
        
        [H]=createGaussianFilter([sigma sigma],[11 11]);
        smoothedImg=filterSeperable(H,img);

        gradMag = calcGradientMag(smoothedImg); 
    end
end